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Abstract. We compute the exact 4-spinon contribution to the zero-temperature 
dynamical structure factor of the spin- 1/2 Heisenberg isotropic antiferromagnet 
in zero magnetic field, directly in the thermodynamic limit. We make use of the 
expressions for matrix elements of local spin operators obtained by Jimbo and 
Miwa using the quantum affine symmetry of the model, and of their adaptation 
to the isotropic case by Abada, Bougourzi and Si-Lakhal (correcting some overall 
factors). The 4-spinon contribution to the first frequency moment sum rule 
at fixed momentum is calculated. This shows, as expected, that most of the 
remaining correlation weight above the known 2-spinon part is carried by 4- 
spinon states. Our results therefore provide an extremely accurate description 
of the exact structure factor. 
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1. Introduction 

It has now been 75 years since Hans Bethe published his seminal paper [I] constructing 
the eigenfunctions of the spin- 1/2 Heisenberg chain [2], 

N 

H = ^S r S j+1 , (1) 

3=1 

giving birth to the Bethe Ansatz and paving the way for the modern theory of 
integrable models of quantum mechanics and field theory [3l [4j [5] . Interest in the 
Heisenberg model has only increased since those early days, partly because of the 
extremely rich mathematical structures now known to be associated to it, but also 
because of its ability to accurately describe a number of real compounds. 

The Bethe Ansatz is first and foremost a method giving access to an integrable 
system's energy levels, allowing the calculation of many equilibrium quantities. For the 
specific case of the Heisenberg model, the ground state energy of the infinite chain was 
computed analytically [6] not long after Bethe's paper, but it was not until the 1960's 
that significant new results were obtained: its excitation spectrum was computed by 
des Cloizeaux and Pearson [7J, and its more general thermodynamic properties were 
obtained shortly afterwards [51 151 HOI HT] . 

Equilibrium quantities are however not sufficient to completely characterize the 
physics of models such as ([1]). Motivated mainly by experimental work, another object 
of fundamental importance has been extensively studied: the dynamical structure 
factor (DSF) 

S ab (k,u;) = Y / e lkl / dte™*<S? +I (t)S*(0)), a,b = x, y, z. (2) 

i=i J ~°° 

This quantity is directly accessible experimentally through inelastic neutron scattering 
[T2l [T3l rm [TB] , and its theoretical calculation opens the door to the interpretation of 
a wealth of experimental data. 

Despite much effort, the dynamics of integrable models remains in general 
inaccessible, and exact calculations based on the Bethe Ansatz can only rarely be 
carried out. One reason is that excitations are usually rather complicated: for the 
Heisenberg model, the Hilbert space can be spanned with spinons [THl H7J [THJ , which 
are nontrivially interacting spin- 1/2 particles whose dispersion relation is given by 

e(p) = ||smp|, pe[-ir,0]. (3) 

The computation of dynamical quantities such as the structure factor however requires 
knowledge of matrix elements of spin operators between (multi-)spinon states, which 
goes much beyond what is accessible with the basic Bethe Ansatz. 

In view of the difficulty of this task, a number of approximate schemes have been 
put forward to offer a qualitative picture of the DSF. One extremely useful construction 
is known as the Miiller ansatz [115] , which is based on exact results for the XY model, 
numerical computations on small chains and known sum rules. Its success lies in its 
extreme simplicity, coupled with rather accurate reproduction of a number of features 
(like the square-root singularity at the lower threshold). It is commonly used in the 
interpretation of experimental data. It drawback is that it is inexact; in particular, 
its functional form at the top of the 2-spinon continuum is not correct. 
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Another important and successful approach relies on mapping the infinite chain 
onto a relativistic quantum field theory [JO] [5T] . Finite size scaling connects the critial 
exponents of the system with its behaviour in a finite volume [HJ [22] > while conformal 
field theory [24, 25 and bosonization allow the calculation of asymptotics of correlation 
functions [26] [27l [28] even at finite temperatures, with known normalizations for the 
first few leading terms in the operator expansion 29, 30, 3 lj . 

As far as methods based on integrability are concerned, it is now possible to 
achieve extremely accurate computations of the DSF over the whole Brillouin zone 
[321 133] for integrable Heisenberg spin chains of any anisotropy in any magnetic 
field, using the ABACUS method [34]. These computations rely on determinant 
representations for matrix elements of local spin operators obtained by solving 
the quantum inverse problem [35[ 136] . and are therefore limited to finite (albeit 
large) chains. We here wish however to concentrate on an altogether different and 
independent take on the problem of calculating the DSF of the Heisenberg model, 
which is not applicable in general but only in one particular (albeit extremely 
important) case. 

The crucial development (as far as the subject of the present paper is concerned) 
came with the recognition that the infinite chain in zero magnetic field displayed a 
quantum affine symmetry, allowing the diagonalization of the Hamiltonian directly 
in the thermodynamic limit [37]. Multi-spinon states then provide a basis for the 
Hilbert space and a resolution of the identity operator, allowing to write the DSF 
as a sum of matrix elements of local spin operators. These matrix elements can be 
computed within this framework through bosonization of the quantum affine algebra 
[38] 139], [40] , a task performed by Jimbo and Miwa [41] . The general representation of 
such correlation functions as the DSF is then obtained in terms of rather complicated 
contour integrals, the number of integrals increasing with the number of spinons in 
the excited state. These are however notoriously hard to evaluate quantitatively. 

One exception is the contribution to the transverse zero-temperature DSF of the 
Heisenberg model coming from 2-spinon intermediate states, the simplest excitations 
that can be constructed above the ground state. In this case, the dynamical constraints 
of conservation of energy and momentum give two 5 functions taking care of the two 
contour integrals involved, allowing for a direct analytical calculation [42] (a similar 
calculation is also possible in the whole gapped antiferromagnetic phase [43]). These 
2-spinon intermediate states were shown to contribute 72.89% of the total structure 
factor intensity [44]. The missing part is then necessarily carried by excited states 
with a higher number of spinons, starting with 4-spinon states. 

The integral representations for matrix elements involving 4-spinon states are 
much more difficult to tackle, since the dynamical constraints are not sufficient to 
take care of all the contour integrals involved (this is also true for the longitudinal 
structure factor of the XX Z model, which was studied in the fsing limit in [45 ). 
The first attempt to tackle more than two spinons was made in [46] . offering a formal 
representation for the exact n-spinon contribution for the zero-temperature DSF of 
the Heisenberg model in zero field. A more thorough treatment of the four-spinon case 
was published shortly afterwards [47], yielding expressions for the DSF in the whole 
gapped antiferromagnetic regime, and their specialization to the Ising and Heisenberg 
limits. These expressions remained however quite complicated, giving the four-spinon 
part at fixed momentum and energy in terms of a double integral of an infinite series. 
Further analytical work [48] yielded little progress, and in fact (as we will show 
below) incorrectly identified the boundaries of the four-spinon continuum. To this 



The 4-spinon dynamical structure factor of the Heisenberg chain 4 

day, nobody has been able to extract curves from these expressions, and previous 
attempts gSJ |3D] have not yielded acceptable results due to the inappropriateness of 
the chosen method and the incorrect continuum used. 

The present paper offers the first reliable computation of the exact 4-spinon 
contribution to the zero-temperature dynamical structure factor of the Heisenberg 
isotropic antiferromagnet in zero magnetic field. After summarizing results for the 
matrix elements which we need for our purposes, the known 2-spinon results are 
repeated as a warmup, after which we present their extension to the case of four 
spinons (correcting some factors in the formulas present in the literature). We finish 
by a discussion of our results, in particular concerning contributions to sum rules. 

2. Exact representation of the dynamical structure factor 

The inevitable first step in the calculation of the DSF © is to insert a resolution 
of the identity in a judiciously chosen basis between the two spin operators. For the 
XXZ model in the thermodynamic limit, a basis for the Hilbert space is provided by 
(multi-) spinon states |£i, £ n )e u ...,e n ;i, which diagonalize the Hamiltonian according 
to 

^l,...,€»>e 1 ,..., en ;i = X;efe)lfl,...,en) 6l ,..., 6 » i < (4) 

where i labels the two equivalent vacuum states |0)j, i = 0, 1 corresponding to the 
two different possible boundary conditions, and = ±1 labels the spin projection of 
the spinons. The coefficients £j are spectral parameters determining the energy and 
momenta of the spinons. The completeness relation reads 

1= EE E i/'n^H ei • ■ : - - • - ^ (5) 

i=0,ln>Oei,..,E„=il ' J j = l * 3 

and can be substituted in ([2]) to yield the decomposition of the DSF into a sum over 
(even numbers of) spinon contributions (here and in the following, we make use of 
spin isotropy and compute S(k,u>) = S zz (k,oj) — S xx = S vy ) 

S{k,u)= ]T S n {k,io). (6) 

n even 

Each term in this decomposition is explicitly written as 

' ra6Zei,...,£„=±l J j=l SJ j=l 

x^Ol^O)!^,...,^,...,^ i;e u ..„e n (6, ...,6,1^(0)10),. (7) 

The matrix elements of local spin operators in the above expression are represented 
exactly within the framework of Jimbo and Miwa [41] . Their adaptation to the 
isotropic Heisenberg antiferromagnet was given in [15] for 2-spinon intermediate states, 
and in [13] for n > 2, although the expressions obtained there are not thoroughly 
simplified. For 4-spinon intermediate states, the matrix elements were studied more 
extensively in [47], whose results will form the basis of the new results we obtain. Let 
us however start by briefly reminding the reader of known results on the much simpler 
case of 2-spinon intermediate states. 
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3. 2-spinon contribution to the structure factor 

It is well-known that 2-spinon intermediate states live within a continuum in k,u) 
defined by satisfying the kinematic constraints of momentum and energy conservation, 

k = -pi-p 2 , oj = e(pi) + e(p 2 ). (8) 

In other words, for a fixed external momentum, there exists an interval in frequency 
given by the conditions 

u > LJ2,i(k) = — | sinfc|, uj < iL>2. u (k) — sin k£[0,2n]. (9) 

The lower boundary is thus given by the des Cloizeaux-Pearson dispersion relation. 
The 2-spinon part of the DSF will be nonvanishing within this continuum, and will by 
construction vanish identically outside of it. This contribution was obtained in [42] 
and is explicitly written as 

S 2 (k,cu) = ~ 6 =e(uj 2iU (k) - u)Q(u - u 2 ,i(k))- (10) 

2 ^i u (k)-^ 

The parameter p is defined as 



coshMM)) = y ^_^ (fc) (ii) 

and the nontrivial part of the DSF is encoded in the fundamental integral function 

m= r^co^)co S y)-i 

7 t cosh<sinh(2t) v ; 

A careful study of this representation of the 2-spinon DSF was carried out in |44| . 
and it is worthwhile to remind the reader of some important facts obtained there. First 
of all, at the lower boundary (for q ^ it), the 2-spinon DSF diverges as a square root 



accompanied by a logarithmic correction, So ~ ; \ /In — ttt (for q = tt, 

the divergence is ~ ^ yjln-^) . Near the upper boundary, on the other hand, the 

2-spinon DSF vanishes in a square-root cusp, S 2 ~ -\Z W2 ."(^) — w - 

Second, sum rules were also studied, most importantly the contribution to the 
total integrated intensity 

•■■*dk rdw i 

o ^io ^' W) = i (13) 
and to the exactly known first frequency moment at fixed momentum, 

Ki{k) = ^ puS(k,u;) = (1 - cosfc)^ (14) 
Jo *™ 

where eo = 1/4 — In 2 is the ground-state energy density [B] . 2-spinon intermediate 
states were shown to carry 72.89% of the total intensity, and 71.30% of the first 
moment sum rule (independently of k). 

More than a quarter of the exact DSF is therefore missing if we restrict ourselves 
to only two spinons. To achieve better saturation of the sum rules, we need to go 
to more complicated intermediate states involving more particles, and we can safely 
expect that out of those, 4-spinon states will be dominant. 
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4. 4-spinon contribution to the structure factor 

Starting from the results of [47j . we write the exact representation for the 4-spinon 
part of the DSF as 

f° 

Si(k, w) = C 4 / dpi...dpi S^ik +} J Pi) S{uj + 22 e{pi)) J({p}) (15) 



where the prefactor is 

3 x 2 9 r(l/4) 8 |yl(i7r/2)| 



e' sinh 2 ft[l + t f]) ^ 
A(z) = exp -/ dt- — — - T — f- . (17) 



with 



t sinh(2£) coshi 

Here again, we restrict to fc G [0, 7r]. Parametrizing the momenta as 

cot = sinh(27rpi) (18) 
the correlation weight is explicitly given by 

J(M) = J({p}) = e-E^WCw) £ | ff;({p}) |2 (19) 

where p^ = p^ — pj. The function /(p) is given by equation p2[) . whereas 5; is given 
by the following expression: 
4 

.9/ = (-iy +1 5>osh(2^-) x 

y, [M^ - W - Q + i Pji ) * r(m - 1 + i^Q 

X m =^-/) n^sinh(^) U.r(m+l + i Pji ) 

where the Heaviside function is here defined as O(n) = for n < and 0(n) = 1 for 
n > 0. 

We have corrected two inaccuracies in [J7]: first, the correct normalization is 
presented here (compare (|20p with formula (5.10) there), and most importantly, 
we have explicitly written that the momentum 8 function fixes k only modulo 
2-7T. This has the crucial consequence that two sectors must be considered when 
solving the dynamical constraints of momentum and energy conservation for 4-spinon 
intermediate states (this was overlooked in (48l [49j |50], leading in particular to an 
incorrect description of the four-spinon continuum) . Bearing in mind that the spinon 
momenta pi are by definition constrained to the interval [— 7r, 0], we define sectors 
and 1 as 

: k+pi +p 2 +p 3 +Pi = 0, 1 : k + 2n + Pl + p 2 + p 3 + Pi = (21) 

with in both cases the energy constraint explicitly written as 

7T 

oj + — (sinpi + sinp2 + sinp3 + sinp4) = 0. (22) 

For higher spinon numbers, more sectors must be similarly added: there are n such 
sectors for states with 2n spinons. 
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A more physical representation of the four-spinon part (|15p of the structure factor 
is obtained by the change of variables {pi} — > {k, uj, K, fi} where 

7T 

K = -pi -p 2 , ft = --(sinpx + sinp 2 )- (23) 

In sector 0, we then have 

7T 

K = k+p 3 +p 4 , ft = uj + -(sinp 3 + sinp 4 ). (24) 

The complete transformation reads 

k fi k n 

Pi = — — + acos — — -, p 2 = — — 



2 oj 2 ,u( K )' y 2 W2,„W 

if-jfc u-Q K-k cj-Sl 
Pa = + acos — , p 4 = — acos — , 25) 

where we restrict to p\ > p 2 and p 3 > p 4 by symmetry. This sector corresponds to 

Ke[o,k}. 

In sector 1, we have 

K = k + 2n+p 3 +p 4 (26) 
instead of the left of (EM)) , yielding the same expressions for pi and p 2 , and 

K-k oj-n K-k Lu-n 

P3 = -7T + — — + aCOS -— -r, Pi = -H+ 7, acos 777 77, (27) 

2 U2, u {K-k) 2 LU2,u(K~k) 

with the same restrictions as above. This sector corresponds to K £ [k, 2tt]. 

The 4-spinon continuum in the k,u> plane is obtained by letting K and Q take 
on all their allowed values. The lower boundary coincides with that of the 2-spinon 
continuum (the des Cloizeaux-Pearson dispersion relation; this is clearly the case from 
a simple physical argument, namely that we are dealing with a massless theory, and 
therefore adding two more spinons of zero momentum in the intermediate state docs 
not shift the energy; in fact, the lower boundary of any Unite higher spinon number 
continuum is identical to the 2-spinon one). The upper boundary, on the other hand, 
extends above the upper boundary of the 2-spinon continuum, and is obtained by 
sharing the momentum (modulo 2n) evenly among the four spinons. Explicitly, we 
therefore have 

W4,i(A) =wa,i(*) = ^| sirxfe), w 4>u (fc) = tt^2 (\ + |cos^. (28) 

A geometrical picture of the 2n spinon continuum is easily obtained by generalization: 
the lower boundary will always be given by the des Cloizeaux-Pearson dispersion 
relation, whereas the upper boundary will be given by the upper boundary of a 2- 
spinon continuum rescaled in size by a factor of n, modulo all its 2tt translations 
(refer to Figure [2] for the simplest case of 4 spinons) . 

For the implementation of the computation of the 4-spinon part of the structure 
factor, it is desirable to describe more carefully the actual integration regions for K 
and £1, which can be precisely defined for given k € [0,7r] and uj € [(^4,i(k),U4, >u (k)]. 
Simple reasoning shows that the K, Q integration regions are obtained by intersecting 
two separate 2-spinon continua, one upright and the other inverted in frequency, and 
shifted with respect to each other by k in momentum and uj in frequency (this is 
illustrated in Figure [1]). The important line crossings are respectively of the u)2. « lines 
(upper boundaries) of the two continua, and of their 0^2, ; lines. The first determine 
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which intervals of K should be included, and the second determine which sub-intervals 
within these should be excluded. 

Regions of K to be included depend on the values of k and u, and are given by 
k 

uj < tt sin - : K G [0, 2tt] 



tt sin \ < uj < 2tt sin | : Ke [if f a , K+] U [K~ b , K+] ( 2! ) I 



where 



* i, - ~ + t ± 2acos- F , Kf b = - ± 2acos— — - . (30) 

2 2tt cos | 2 2tt sm | 

The lower boundaries of the intersecting continua define excluded regions of K as 
| sin k < uj < tt sin - : X g [JC£, if +] U [if 2 ~ + tt, if + + tt], 

|sinfc< W <7rcos| : if £ [if^, if+ ] U [if^ + tt, K+ + tt] (31) 

where 

= TJ + 77 ± acos ~ k . = J iacos^^. (32) 

2 2 7r cos | 2 7r sm f 

We define the if integration domain T>k as the sum of (possibly disconnected) regions 
fulfilling the above constraints. 

For fixed k, uj, if fulfilling the above constraints, the value of f2 is restricted to a 
finite interval: 

Qi(k, uj,K) <fl< fl u (k, uj, if), (33) 
with the limits being explicitly given by 

f tt k — if 

Q;(fc, uj, if) = Max ( — | sin if |, uj — tt sin | — - — 

n u (k,uj,K) = Min ^7rsiny,w- ||sin(fc- K)\j . (34) 

The leftover two-dimensional integral for the 4-spinon part of the DSF is therefore 
over a region with nontrivial geometry, depending on the particular values of k and 
uj. Within the 4-spinon continuum, we can identify six sectors (illustrated in Figure 
[2]), each of which leads to a different sort of integration domain in the if, plane 
(examples of which are illustrated in Figured]). By symmetry, we only need to consider 
k G [0,tt]. 

In terms of the new variables k,uj, K,£l, the kinematic restrictions are trivially 
implemented, and the four-spinon contribution to the structure factor can be written 
as a two-dimensional integral in if and f2 over the regions defined above: 

Si(k,oj) = c 4 f dK [ nuik ' a ' K dn J(k,uj,K,n) 

' ' - la lik ,u,K) { [ujI jk) - w] [aj 3 u{k - K )-{uj- ny] 

where we have written J implicitly as a function of the new variables. 

For the evaluation of this expression, we use a specially adapted Romberg- 
like integration method. The integrand typically has integrable divergences at the 
boundaries of the integration regions, which are absorbed by various appropriate 
changes of variables before coding. While technically feasible as we demonstrate 
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Figure 1. Integral regions in the K, O plane. These are formed by the intersection 
of two 2-spinon continua, one of them inverted, and shifted with respect to one 
another by k and ui. The left-hand example is a simple case where only one 
connected domain is obtained. On the right is a more complicated example with 
three different integration domains (whose K borders are pointed out above; the 
domains delimited by a and c have a very small but finite area). 




Figure 2. For values of k 
and ui within the 4-spinon 
continuum (shaded), the 
K and Q integration re- 
gions take different forms. 
Six sectors are obtained, 
as labeled here (sector 5 is 
barely visible, bordered by 
sectors 2,3,6). Examples 
of K, Q integration regions 
for each of these sectors are 
given in Figure [4] 






A 



Figure 3. Integration domains in the K, Q plane (shaded) for the six sectors 
given in Figuref2](l — 3 top left to right, 4 — 6 bottom). The K = 0, Q = origin 
always lies at the left foot of the upright 2-spinon continuum. 
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Figure 4. Plots of the 2- and 4-spinon parts of the dynamical structure factor 
at momenta k = it /A, n/2, 3n/4 and n. The horizontal axis is frequency, and we 
focus on the region of the 2-spinon continuum for ease of comparison. 



here, it remains a challenge to obtain very good precision, and we therefore here 
first concentrate on results for fixed momentum. 

The main results are plotted in Figure 2J for four representative values of k 
(other values give similar-looking plots). Over the 2-spinon continuum, the 4-spinon 
contribution is of the same order as the 2-spinon one (i.e. about a third of it). Between 
the upper boundary u>2,u of the 2-spinon continuum and the upper boundary io^ u of 
the 4-spinon continuum, the 4-spinon part of the structure factor is finite but very 
small. Figure 5 displays the shape of the DSF in the vicinity of the 2-spinon upper 
boundary for two representative values of momentum (the same sort of behaviour is 
observed for all momentum values we checked). A few things are worth pointing out 
here. First, at the lower boundary, the 4-spinon part diverges similarly to the 2-spinon 
one. Second, the 4-spinon part is finite and smooth around the upper boundary W2, u 
of the 2-spinon continuum. The full DSF therefore still has a square-root singularity 
around this point, yielding a picture consistent with that put forward in [51j (see 
also the related discussion in [52]). At higher frequencies, the DSF decays extremely 
rapidly, as illustrated in Figure [6] 

To assess the quality of our results at fixed momentum, we compute the first 
frequency moment sum rule in Table 1. 4-spinon intermediate states clearly carry (as 
expected) the majority of the missing correlation weight after 2-spinon contributions 
have been taken into account, i.e. around 27% ±1% of this sum rule for all values 
of momentum which were studied. We have not yet achieved sufficient accuracy 
for the 4-spinon part of the total integrated intensity, but we can expect again a 



The 4-spinon dynamical structure factor of the Heisenberg chain 



11 




Figure 5. Zooms on the upper boundary of the 2-spinon continuum, here for 
k = 7r/4 and ir/2 (other values of momentum give very similar results). The 
2-spinon part vanishes in a square root cusp at u>2, u , whereas the 4-spinon part 
is finite. There is thus an infinite slope in the full dynamical structure factor at 



-km It/4 

-k = % 



1 2 3 4 5 6 
CO 



Figure 6. The 4-spinon 
part (in logarithmic scale) 
of the DSF as a function of 
frequency for three repre- 
sentative values of momen- 
tum. The upper boundary 
of the 2-spinon continuum 
is visible as the shoulder 
in the curve, with the 4- 
spinon part of the DSF de- 
creasing rapidly above this 
threshold. 



k 


tt/8 


tt/4 


37r/8 


tt/2 


5tt/8 


3tt/4 


7tt/8 


7T 


K x (k) (4 s P %) 


27.2 


27.3 


27.3 


27.1 


26.5 


26.2 


25.9 


26.6 



Table 1. First frequency moment sum rule (see equation JH}) percentages 
coming from the 4-spinon part of the dynamical structure factor. The 2-spinon 
part always contributes 71.30%. 4-spinon states clearly carry the majority of the 
leftover correlation weight, as is naturally expected. The accuracy of these results 
is estimated to be around 1% of K\{k). 



contribution of the order of 27% (since the general shape of the 4-spinon part resembles 
that of the 2-spinon part, and since the relative 2-spinon part of both sum rules 
is almost equal). There thus remains only about 2% missing, which are naturally 
ascribed to higher spinon numbers. Although these are also in principle accessible 
using an extension of the present method, the (for 6 spinons, quadruple) integrations 
needed probably prohibit accurate evaluation without first achieving significant further 
analytical advances. However, this missing part is now rather small, meaning that our 
results provide an approximation of the exact zero-temperature DSF of the Heisenberg 
model in zero field which is the most accurate available at the moment. 
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Conclusion 

In conclusion, we have calculated the 4-spinon contribution to the zero-temperature 
dynamical structure factor of the Heisenberg model in zero magnetic field, starting 
from its exact integral representation in the thermodynamic limit. As proven by 
sum rules, this contribution carries most of the correlation weight left over after 2- 
spinon intermediate states have been taken into accout. The results obtained therefore 
provide a very close description of the exact correlator. In future publications, we will 
provide a thorough analysis of other available sum rules, and extend the results to the 
gapped antiferromagnetic case. 
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